#######################################################################################
### R file to replicate the figures in the WWW JoP project
###
### Created: 4-9-15
### Modified: 9-5-19
###
#######################################################################################

library(foreign)
library(ggplot2)
library(lattice)
library(fields)
library(splines)
library(MASS)
library(gridExtra)

slx = read.dta("Figures/SLX.dta", convert.underscore=TRUE)

#######################################################################################
### Figure 2: Recovery of the True Total Effects of X1 in the SAR Model across Strength of Spatial Autocorrelation and Spatial Heterogeneity
#######################################################################################
ls <- 1.5
bs <- 30

id1 <- ggplot(subset(slx, id == 1), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.2), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = -0.8"))) + xlab("") + ylab("") + theme_minimal()
id2 <- ggplot(subset(slx, id == 2), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.2), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = 0.01"))) + xlab("") + ylab("") + theme_minimal()
id3 <- ggplot(subset(slx, id == 3), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.2), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()
id4 <- ggplot(subset(slx, id == 4), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.8), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = -0.2"))) + xlab("") + ylab("") + theme_minimal()
id5 <- ggplot(subset(slx, id == 5), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.8), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = 0.01"))) + xlab("") + ylab("") + theme_minimal()
id6 <- ggplot(subset(slx, id == 6), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.8), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()
id7 <- ggplot(subset(slx, id == 7), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 1.4), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = -0.8"))) + xlab("") + ylab("") + theme_minimal()
id8 <- ggplot(subset(slx, id == 8), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 1.4), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = 0.01"))) + xlab("") + ylab("") + theme_minimal()
id9 <- ggplot(subset(slx, id == 9), aes(teffect.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 1.4), size = ls) + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()

# Combine the plots onto one figure
grid.newpage()
pushViewport(viewport(layout = grid.layout(3,3)))

vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(id1, vp = vplayout(1, 1))
print(id2, vp = vplayout(1, 2))
print(id3, vp = vplayout(1, 3))
print(id7, vp = vplayout(2, 3))
print(id8, vp = vplayout(2, 2))
print(id9, vp = vplayout(2, 1))
print(id4, vp = vplayout(3, 1))
print(id5, vp = vplayout(3, 2))
grid.text(expression(paste("Average Total Effects, ", x[1])), x=.5, y=0.075, vp = vplayout(3, 2))
print(id6, vp = vplayout(3, 3))



#######################################################################################
### Figure 3: Recovery of the True Indirect Effects of X2 in the SAR Model across Strength of Spatial Autocorrelation and Spatial Heterogeneity
#######################################################################################
ls <- 1.5
bs <- 30

id1 <- ggplot(subset(slx, id == 1), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = -0.8), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = -0.8"))) + xlab("") + ylab("") + theme_minimal()
id2 <- ggplot(subset(slx, id == 2), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.01), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = 0.01"))) + xlab("") + ylab("") + theme_minimal()
id3 <- ggplot(subset(slx, id == 3), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.4), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()
id4 <- ggplot(subset(slx, id == 4), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = -0.2), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = -0.2"))) + xlab("") + ylab("") + theme_minimal()
id5 <- ggplot(subset(slx, id == 5), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.01), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = 0.01"))) + xlab("") + ylab("") + theme_minimal()
id6 <- ggplot(subset(slx, id == 6), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.4), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()
id7 <- ggplot(subset(slx, id == 7), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = -0.8), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = -0.8"))) + xlab("") + ylab("") + theme_minimal()
id8 <- ggplot(subset(slx, id == 8), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.01), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = 0.01"))) + xlab("") + ylab("") + theme_minimal()
id9 <- ggplot(subset(slx, id == 9), aes(ieffect.x2.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0.4), size = ls) + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(ieffect.x2.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()

# Combine the plots onto one figure
grid.newpage()
pushViewport(viewport(layout = grid.layout(3,3)))

vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(id1, vp = vplayout(1, 1))
print(id2, vp = vplayout(1, 2))
print(id3, vp = vplayout(1, 3))
print(id7, vp = vplayout(2, 3))
print(id8, vp = vplayout(2, 2))
print(id9, vp = vplayout(2, 1))
print(id4, vp = vplayout(3, 1))
print(id5, vp = vplayout(3, 2))
grid.text(expression(paste("Average Indirect Effects, ", x[2])), x=.5, y=0.075, vp = vplayout(3, 2))
print(id6, vp = vplayout(3, 3))


#######################################################################################
### Figure 4: SAR Model Indicates the Presence of Higher-Order Indirect Effects for x1 When They Are Not There
#######################################################################################
ls <- 1.5
bs <- 30

id1.o2 <- ggplot(subset(slx, id == 1), aes(teffect.o2.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0), size = ls) + geom_vline(aes(xintercept = quantile(teffect.o2.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.o2.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = -0.8"))) + xlab("") + ylab("") + theme_minimal()
id9.o2 <- ggplot(subset(slx, id == 9), aes(teffect.o2.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0), size = ls) + geom_vline(aes(xintercept = quantile(teffect.o2.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.o2.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()
id4.o2 <- ggplot(subset(slx, id == 4), aes(teffect.o2.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0), size = ls) + geom_vline(aes(xintercept = quantile(teffect.o2.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.o2.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = -0.2"))) + xlab("") + ylab("") + theme_minimal()

id1.o3 <- ggplot(subset(slx, id == 1), aes(teffect.o3.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0), size = ls) + geom_vline(aes(xintercept = quantile(teffect.o3.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.o3.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.8; ", theta[2], " = -0.8"))) + xlab("") + ylab("") + theme_minimal()
id9.o3 <- ggplot(subset(slx, id == 9), aes(teffect.o3.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0), size = ls) + geom_vline(aes(xintercept = quantile(teffect.o3.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.o3.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = 0.4; ", theta[2], " = 0.4"))) + xlab("") + ylab("") + theme_minimal()
id4.o3 <- ggplot(subset(slx, id == 4), aes(teffect.o3.x1.m4)) + geom_histogram(bins = bs, fill=I("grey50")) + geom_vline(aes(xintercept = 0), size = ls) + geom_vline(aes(xintercept = quantile(teffect.o3.x1.m4, c(.025))), linetype="dashed") + geom_vline(aes(xintercept = quantile(teffect.o3.x1.m4, c(.975))), linetype="dashed") + ggtitle(expression(paste(theta[1], " = -0.2; ", theta[2], " = -0.2"))) + xlab("") + ylab("") + theme_minimal()

# Combine the plots onto one figure
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,3)))

vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(id4.o2, vp = vplayout(1, 1))
print(id9.o2, vp = vplayout(1, 2))
print(id1.o2, vp = vplayout(1, 3))

print(id4.o3, vp = vplayout(2, 1))
print(id9.o3, vp = vplayout(2, 2))
print(id1.o3, vp = vplayout(2, 3))